iHMP

library(data.table)
library(magrittr)
library(ggplot2)
theme_set(theme_minimal())

metadata <- fread("ihmp/data/ibdmdb_healthy.csv")
sra_data <- fread("ihmp/data/ibdmd_healthy_runtable.tsv")[, .(Run, `Library Name`)]
sra_data[, "External ID" := tstrsplit(`Library Name`, "_MGX")[[1]]]
metadata <- metadata[sra_data, on="External ID", nomatch=0]
metadata

man <- fread("db/data/manifest.csv")[, .(db, rank, matched_taxid, seqlength)]

foods <- fread("ihmp/data/food_abundance.csv")
foods <- man[foods, on="matched_taxid"]
foods[, "tpm" := 1.0e6 * reads / as.double(seqlength)]
foods[is.na(tpm), tpm := 0.0]

contents <- fread("ihmp/data/food_content.csv")
contents

tests <- list(
  vegetables = list(group="Vegetables (salad, tomatoes, onions, greens, carrots, peppers, green beans, etc)", ex=expression(food_group == "Vegetables")),
  fruits = list(group="Fruits (no juice) (Apples, raisins, bananas, oranges, strawberries, blueberries", ex=expression(food_group == "Fruits")),
  beans = list(group="Beans (tofu, soy, soy burgers, lentils, Mexican beans, lima beans etc)", ex=expression(food_subgroup == "Beans")),
  `white meat` = list(group="White meat (chicken, turkey, etc.)", ex=expression(food_subgroup == "Poultry")),
  shellfish = list(group="Shellfish (shrimp, lobster, scallops, etc.)", ex=expression(food_subgroup %chin% c("Mollusks", "Crustaceans"))),
  fish = list(group="Fish (fish nuggets, breaded fish, fish cakes, salmon, tuna, etc.)", ex=expression(food_subgroup == "Fishes")),
  `red meat` = list(group="Red meat (beef, hamburger, pork, lamb)", ex=expression(food_subgroup %chin% c("Swine", "Bovines", "Caprae")))
)

library(patchwork)

tables <- list()

freqs <- c(
  `No, I did not consume these products in the last 7 days` = 0,
  `Within the past 4 to 7 days` = 1/5.5,
  `Within the past 2 to 3 days` = 1/2.5,
  `Yesterday, 1 to 2 times` = 1.5,                                
  `Yesterday, 3 or more times` = 3,
  `NA` = 0
)

results <- list()
plots <- list()
merged <- metadata[foods, on=c(`Run` = "sample_id")]
for (i in 1:length(tests)) {
  name <- names(tests)[i]
  vals <- tests[[i]]
  dt <- merged[, .(abundance=sum(reads[eval(vals$ex)]), group=.SD[[vals$group]][1], total_reads=total_reads[1]), by="External ID"]
  ns <- dt[, .N, by="group"] |> setkey(group)
  dt <- dt[group %chin% names(freqs) & ns[group, N] >= 10]
  dt[, "frequency" := freqs[group]]
  dt[, group := factor(group, levels=names(freqs)[names(freqs) %in% group])]
  dt[, "relative" := (abundance + 1) / (total_reads)]
  dt[, "item" := name]
  plots[[name]] <- ggplot(dt) + aes(x=group, y=log10(relative)) + 
    geom_jitter(width=0.3) + 
    geom_boxplot(width=0.1) + 
    stat_smooth(aes(group=1), method="lm") + 
    labs(title=name, x="consumption frequency [servings/day]", y="abundance")
  print(name)
  abmod <- lm(log10(relative) ~ frequency, data=dt)
  premod <- glm((abundance > 0) ~ frequency, data=dt, family=binomial(link="logit"))
  tables[[name]] <- dt
  results[[name]] <- dt[, .(n=.N, relative=mean(log10(relative*1e6)), sd=sd(log10(relative*1e6)), detected=sum(abundance > 0), item=name), by="group"]
  results[[name]][, "p_abundance" := summary(abmod)$coefficients[2, 4]]
  results[[name]][, "p_prevalence" := summary(premod)$coefficients[2, 4]]
}

[1] "vegetables"
[1] "fruits"
[1] "beans"
[1] "white meat"
[1] "shellfish"
[1] "fish"
[1] "red meat"

r <- rbindlist(results)
ta <- rbindlist(tables)
ggplot(r) +
  aes(x=group, y=relative) +
  geom_jitter(data=ta, aes(y=log10(relative*1e6)), stroke=0, size=0.5, width=0.3) +
  geom_pointrange(aes(ymin=relative-sd, ymax=relative+sd), shape=21, fill="white") +
  geom_text(data=unique(r, by="item"), aes(x=0, y=Inf, label=sprintf("p[lm]=%.2g\np[logit]=%.2g", p_abundance, p_prevalence)), size=3, vjust=1, hjust=0, nudge_y=0.5) +
  facet_wrap(~ item, ncol=2) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  labs(x="", y=expression(paste("log relative abundance [",log[10],RPM, "]"))) +
  guides(color=F)

ggsave("figures/ihmp_foods.pdf", width=2.5, height=9)

ggplot(foods[seqlength>1e6], aes(x=seqlength/1e9, y=relative_abundance)) + geom_point() + stat_smooth(method="lm", color="tomato") +
  labs(x="haploid genome size [Gbps]", y="relative food abundance") + scale_x_log10() + scale_y_log10()

foods[seqlength>1e6, cor.test(reads, seqlength)]

    Pearson's product-moment correlation

data:  reads and seqlength
t = 1.9748, df = 2209, p-value = 0.04841
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.000294607 0.083521188
sample estimates:
       cor 
0.04198072 

ggsave("figures/genome_size.pdf", width=5, height=3)

Diet distances

species <- fread("ihmp/data/S_counts.csv")[grepl("Bacteria|Archaea", lineage)]
microbiome_mat <- dcast(species, sample_id ~ lineage, value.var = "new_est_reads", fill=0, fun.aggregate = sum)
sids <- microbiome_mat[, tstrsplit(sample_id, "S_")[[2]]]
microbiome_mat <- as.matrix(microbiome_mat[, "sample_id" := NULL])
rownames(microbiome_mat) <- sids
good <- colMeans(microbiome_mat) > 10
microbiome_mat <- microbiome_mat[, good]
# micro_rare <- phyloseq(otu_table(microbiome_mat, taxa_are_rows = F)) |> rarefy_even_depth(100000) |> otu_table()
microbiome_relative <- microbiome_mat / rowSums(microbiome_mat)

diet <- metadata[, names(metadata)[72:92][-11], with=F]
diet <- apply(diet, 2, function(x) ifelse(x == "", 0, freqs[x]))
rownames(diet) <- metadata[["Run"]]
diet <- diet[rowSums(diet) > 0, ]
diet <- diet[, colMeans(diet >0) > 0.1]
dim(diet)

[1] 362  19

foodab <- foods
foodab[, "id" := paste(matched_taxid, species, wikipedia_id, sep="|")]
food_mat <- dcast(foodab, sample_id ~ id, value.var = "reads", fill=0, fun.aggregate = sum)
sids <- food_mat[, sample_id]
food_mat <- as.matrix(food_mat[, sample_id := NULL])
rownames(food_mat) <- sids
food_relative <- food_mat[rowSums(food_mat) > 0, ]
food_relative <- food_relative / rowSums(food_relative)
food_relative <- food_relative[rowSums(food_relative) > 0, ]
food_relative <- food_relative[, colMeans(food_relative >0) > 0.1]
dim(food_relative)

[1] 359  16

Beta diversity tests

Let’s write a little function that runs the test for microbiome <-> other comparison and plots and returns the results.

library(vegan)

micro_mantel <- function(first, other, firstname, othername) {
  sids <- intersect(rownames(first), rownames(other))
  first_dist <- vegdist(first[sids, ], "bray")
  other_dist <- vegdist(other[sids, ], "bray")
  
  test <- mantel(first_dist, other_dist, method="pearson")
  results <- data.table(first=as.numeric(first_dist), other=as.numeric(other_dist))
  pl <- ggplot(results) +
    aes(x=first, y=other) +
    geom_point(size=1, alpha=0.25, stroke=0) +
    stat_smooth(method="lm") +
    labs(x=paste(firstname, "distance [Bray]"), y=paste(othername, "distance [Bray]"))
  print(pl)
  return(data.table(x=firstname, y=othername, r=test$statistic, p=test$signif, perms=tests$permutations))
}

Now we run it for the measures.

mantel_tests <- list(
  micro_mantel(microbiome_relative, diet, "microbiome", "FFQ"),
  micro_mantel(microbiome_relative, food_relative, "microbiome", "MEDI"),
  micro_mantel(diet, food_relative, "FFQ", "MEDI")
) |> rbindlist()

mantel_tests

Figure data

fwrite(ta, "figure_data/iHMP_points.csv")
fwrite(r, "figure_data/iHMP_summaries.csv")
fwrite(mantel_tests, "figure_data/iHMP_mantel.csv")

fwrite(foods[seqlength>1e6], "figure_data/ihmp_genome_size.csv")

Analysis for hibiscus reads

library(patchwork)

tables <- list()

freqs <- c(
  `No, I did not consume these products in the last 7 days` = 0,
  `Within the past 4 to 7 days` = 1/5.5,
  `Within the past 2 to 3 days` = 1/2.5,
  `Yesterday, 1 to 2 times` = 1.5,                                
  `Yesterday, 3 or more times` = 3,
  `NA` = 0
)
tests <- names(metadata)[72:92][-11]

results <- list()
plots <- list()
merged <- metadata[foods, on=c(`Run` = "sample_id")]
for (i in 1:length(tests)) {
  name <- tests[i]
  dt <- merged[, .(abundance=sum(reads[genus=="Hibiscus"]), group=.SD[[name]][1], total_reads=total_reads[1]), by="External ID"]
  ns <- dt[, .N, by="group"] |> setkey(group)
  dt <- dt[group %chin% names(freqs) & ns[group, N] >= 10]
  dt[, "frequency" := freqs[group]]
  dt[, group := factor(group, levels=names(freqs)[names(freqs) %in% group])]
  dt[, "relative" := (abundance + 1) / (total_reads)]
  dt[, "item" := name]
  if (dt[, uniqueN(group) == 1])  {
    next
  }
  plots[[name]] <- ggplot(dt) + aes(x=group, y=log10(relative)) + 
    geom_jitter(width=0.3) + 
    geom_boxplot(width=0.1) + 
    stat_smooth(aes(group=1), method="lm") + 
    labs(title=name, x="consumption frequency [servings/day]", y="abundance")
  print(name)
  abmod <- lm(log10(relative) ~ frequency, data=dt)
  premod <- glm((abundance > 0) ~ frequency, data=dt, family=binomial(link="logit"))
  tables[[name]] <- dt
  results[[name]] <- dt[, .(n=.N, relative=mean(log10(relative*1e6)), sd=sd(log10(relative*1e6)), detected=sum(abundance > 0), item=name), by="group"]
  results[[name]][, "p_abundance" := summary(abmod)$coefficients[2, 4]]
  results[[name]][, "p_prevalence" := summary(premod)$coefficients[2, 4]]
}

[1] "Soft drinks, tea or coffee with sugar (corn syrup, maple syrup, cane sugar, etc)"
[1] "Diet soft drinks, tea or coffee with sugar (Stevia, Equal, Splenda etc)"
[1] "Fruit juice (orange, apple, cranberry, prune etc.)"
[1] "Water"
[1] "Alcohol (beer, brandy, spirits, hard liquor, wine, aperitif, etc.)"
[1] "Yogurt or other foods containing active bacterial cultures (kefir, sauerkraut)"
[1] "Dairy (milk, cream, ice cream, cheese, cream cheese)"
[1] "Fruits (no juice) (Apples, raisins, bananas, oranges, strawberries, blueberries"
[1] "Vegetables (salad, tomatoes, onions, greens, carrots, peppers, green beans, etc)"
[1] "Beans (tofu, soy, soy burgers, lentils, Mexican beans, lima beans etc)"
[1] "Whole grains (wheat, oats, brown rice, rye, quinoa, wheat bread, wheat pasta)"
[1] "Starch (white rice, bread, pizza, potatoes, yams, cereals, pancakes, etc.)"
[1] "Eggs"
[1] "Processed meat (other red or white meat such as lunch meat, ham, salami, bologna"
[1] "Red meat (beef, hamburger, pork, lamb)"
[1] "White meat (chicken, turkey, etc.)"
[1] "Shellfish (shrimp, lobster, scallops, etc.)"
[1] "Fish (fish nuggets, breaded fish, fish cakes, salmon, tuna, etc.)"
[1] "Sweets (pies, jam, chocolate, cake, cookies, etc.)"

r <- rbindlist(results)
ta <- rbindlist(tables)
sig <- r[p_abundance < 0.05/length(results)]
ggplot(sig) +
  aes(x=group, y=relative) +
  geom_jitter(data=ta[item %in% sig$item], 
              aes(y=log10(relative*1e6)), stroke=0, size=0.5, width=0.3) +
  geom_pointrange(aes(ymin=relative-sd, ymax=relative+sd), shape=21, fill="white") +
  geom_text(data=unique(sig, by="item"), aes(x=0, y=Inf, label=sprintf("p[lm]=%.2g\np[logit]=%.2g", p_abundance, p_prevalence)), size=3, vjust=1, hjust=0, nudge_y=0.5) +
  facet_wrap(~ item, nrow=1, labeller = labeller(item = label_wrap_gen(40))) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  labs(x="", y=expression(paste("Hibiscus abundance [",log[10],RPM, "]"))) +
  guides(color=F)

ggsave("figures/hibiscus.pdf", width=8, height=6)

Figure Data

openxlsx::write.xlsx(
  list(`Fig. ED6 (individual)` = ta, `Fig. ED6 (summaries)` = r),
  "figure_data/figED6.xlsx"
)

fwrite(metadata, "figure_data/ihmp_metadata.csv")